辅助工具使用教程 ************************************ .. note:: cli 是 dspawpy （一个基于 :bdg-info:`Python` :bdg-info-line:`≥ 3.8` 的库）的命令行交互程序，提供了十几个选项，用于快速完成常见任务。对于更复杂的任务，可参考本章第二节开始的脚本示例。 :ref:`安装 ` dspawpy后，命令行输入 ``cli`` 即可使用： .. code-block:: text ********这是dspawpy命令行交互小工具，预祝您使用愉快******** ( ) _| | ___ _ _ _ _ _ _ _ _ _ _ _ /'_` |/',__)( '_`\ /'_` )( ) ( ) ( )( '_`\ ( ) ( ) ( (_| |\__, \| (_) )( (_| || \_/ \_/ || (_) )| (_) | \__,_)(____/| ,__/'`\__,_) \___x___/ | ,__/ \__, | | | | | ( )_| | (_) (_) \___/ 1.3.0: /home/zzl_ubuntu22/Docs/dspaw-manual-cn/dspawpy_proj/dspawpy 1: update更新 2: structure结构转化 3: volumetricData数据处理 4: band能带数据处理 5: dos态密度数据处理 6: bandDos能带和态密度共同显示 7: optical光学性质数据处理 8: neb过渡态计算数据处理 9: phonon声子计算数据处理 10: aimd分子动力学模拟数据处理 11: Polarization铁电极化数据处理 12: ZPE零点振动能数据处理 13: TS的热校正能 --> 输入数字后回车选择功能： .. _install instruction: 安装与更新 ============================================================================== 1. 在HZW机器上，已提前安装 dspawpy，使用以下命令激活虚拟环境即可使用： .. code-block:: shell source /data/hzwtech/profile/dspawpy.env 2. 在其他机器上，请自行安装 dspawpy（下列两种方式二选一）： - 使用 conda 推荐使用 `conda `_ 它在管理 python 包时会全局考虑依赖关系，出现底层依赖库错误的概率较低 .. code-block:: shell conda install dspawpy -c conda-forge - 或使用 pip .. dropdown:: :octicon:`code-review;1.5em;sd-text-light` pip :animate: fade-in-slide-down :color: info - pip 是 python 自带的包管理器 - Linux和Mac一般自带python和pip。 - Windows平台打开 microsoft store 搜索 python 安装即可 .. figure:: ../imgs/python_ms_store.png :height: 300px :align: center 然后打开 cmd 或者 powershell 即可使用 pip。 .. code-block:: shell pip install dspawpy 更新 dspawpy --------------------------------------------------- 如果 dspawpy 是通过 conda 安装的，使用以下命令更新： .. code-block:: shell conda update dspawpy 如果 dspawpy 是通过 pip 安装的，那么： .. code-block:: shell pip install dspawpy -U # -U 表示安装最新版 .. note:: 如果 pip 使用了国内的镜像网站，可能由于镜像网站尚未同步最新版 `dspawpy `_ |dspawpy_version| 而导致无法顺利升级，需指定从pypi官网下载安装： .. code-block:: shell pip install dspawpy -i https://pypi.org/simple --user -U # -i 指定下载地址，--user 表示仅为当前用户安装，-U 表示安装最新版 如果运行时出现 **dspawpy 相关** 错误信息，请先检查是否已正确导入最新版本的 dspawpy 并检查安装路径： .. code-block:: shell $ python3 # 或 python >>> import dspawpy >>> dspawpy.__version__ # 将输出版本号 >>> dspawpy.__file__ # 将输出安装路径 如果当前用户目录下存在 ``.local/lib`` 文件夹，Python 将优先从中寻找并导入包，这可能导致包的版本错误（即使已进入conda虚拟环境）！ 请尝试在运行时增加 -s 选项避免导入 .local/lib 中的包： ``python -s your-script.py`` .. _结构转化: structure结构转化 ============================================================================== 读取结构信息使用 ``read`` 函数；将结构信息写入文件，使用 ``write`` 函数；快速结构转化，使用 ``convert`` 函数： .. dropdown:: :octicon:`code-review;1.5em;sd-text-light` API: read(), write(), convert() :animate: fade-in-slide-down :color: info .. automodule:: dspawpy.io.structure :members: read, write, convert :noindex: 可参考 :download:`2conversion.py <../dspawpy_proj/UserScripts/2conversion.py>` 脚本完成转化： .. literalinclude:: ../dspawpy_proj/UserScripts/2conversion.py :language: python :linenos: convert 函数的几个关键参数设置规则见下表： .. list-table:: **dspawpy 支持读写的结构文件列表** :widths: 36 12 14 14 38 :header-rows: 1 * - infile/outfile（输入/输出文件名模糊匹配） - infmt（输入文件格式） - outfmt（输出文件格式） - 兼容性 - 说明 * - `*.h5` - h5 - h5 - 读+ - DS-PAW计算完成后保存的hdf5类型文件 * - `*.json` - json - json - 读+ 写 - DS-PAW计算完成后保存的json类型文件 * - `*.pdb` - pdb - pdb - 读+ 写+ - Protein Data Bank * - `*.as` - as - as - 读 写 - DS-PAW记载原子坐标等信息的结构文件 * - `*.hzw` - hzw - hzw - 读 写 - DeviceStudio默认的结构文件 * - `*.xyz` - 无 - xyz - 写+ - extended-xyz类型轨迹文件 * - `*.dump` - 无 - dump - 写+ - lammps-dump类型轨迹文件 * - `*.cif*/*.mcif*` - 无 - cif/mcif - 读 写 - Crystallographic Information File * - `*.vasp/*POSCAR*/*CONTCAR*` - 无 - poscar - 读 写 - VASP结构文件 * - `*.cssr*` - 无 - cssr - 读 写 - Crystal Structure Standard Representation * - `*.xsf*` - 无 - xsf - 读 写 - eXtended Structural Format * - `*rndstr.in*/*lat.in*/*bestsqs*` - 无 - mcsqs - 读 写 - Monte Carlo Special Quasirandom Structure * - `*.yaml/*.yml` - 无 - yaml - 读 写 - YAML Ain't Markup Language * - `inp*.xml/*.in*/inp_*` - 无 - fleur-inpgen - 读 写 - FLEUR 结构文件，须安装 pymatgen-io-fleur * - `*prismatic*` - 无 - prismatic - 写 - 用于STEM模拟的一种文件格式 * - `*.res` - 无 - res - 写 - ShelX res 结构文件 .. note:: - 读取时，infmt参数优先级高于文件名模糊匹配规则；例如，指定infmt='h5'后，文件名可以是任意值，甚至可以是 a.json - 写入时，outfmt参数优先级高于文件名模糊匹配规则；例如，指定outfmt='h5'后，文件名可以是任意值，甚至可以是 a.json - “infile/outfile”列 ``*`` 号表示任意字符，“兼容性”列 ``+`` 号表示支持处理多个结构，否则只能处理单个结构： .. warning:: - 将结构信息写入文件时，所有列出的格式(outfmt)都支持，而文件名模糊匹配规则不一定完全适用 - `*.cif*/*.mcif*` 及以下的匹配规则已简化，更多参数细节，请查阅 pymatgen 的 `from_file() `_ 和 `to() `_ 两个函数的源代码。 - 将结构信息写为 `json` 格式时，仅支持可视化 NEB 链任务，详见 :ref:`观察NEB链` 节相关说明 - DS-PAW 输出的 neb.h5、phonon.h5、phonon.json、neb.json以及wannier.json 暂时无法读取结构信息。 .. _ 电荷密度数据处理: volumetricData数据处理 ============================================================================== volumetricData可视化 ----------------------------------------------------- - 参考 :download:`3vis_vol.py <../dspawpy_proj/UserScripts/3vis_vol.py>` ： .. literalinclude:: ../dspawpy_proj/UserScripts/3vis_vol.py :language: python :linenos: 将转换后的文件 :guilabel:`DS-PAW_rho.cube` 拖入 **VESTA** 软件中显示： .. figure:: ../imgs/1-3.png :height: 300px :align: center .. centered:: 晶体硅单胞的电荷密度分布示意图 差分volumetricData可视化 ----------------------------------------------------- - 参考 :download:`3dvol.py <../dspawpy_proj/UserScripts/3dvol.py>` ： .. literalinclude:: ../dspawpy_proj/UserScripts/3dvol.py :language: python :linenos: 上述脚本支持处理多元体系的电荷密度差，示例以二元体系为例，得到了 AB.h5 与 A.h5 和 B.h5 之间的电荷密度差文件 :guilabel:`delta_rho.cube` ，该文件可直接使用 **VESTA** 打开。 volumetricData面平均 ----------------------------------------------------- - 参考 :download:`3planar_ave.py <../dspawpy_proj/UserScripts/3planar_ave.py>` ： .. literalinclude:: ../dspawpy_proj/UserScripts/3planar_ave.py :language: python :linenos: 处理 :doc:`/appliedcases` 3.3 小节所得静电势文件，可得真空方向势函数曲线如下所示： .. figure:: ../imgs/us/3pot_ave.png :height: 450px :align: center .. centered:: AuAl势函数示意图 .. dropdown:: :octicon:`code-review;1.5em;sd-text-light` API: write_VESTA(), write_delta_rho_vesta(), average_along_axis() :animate: fade-in-slide-down :color: info - ``write_VESTA`` 函数负责处理volumetricData可视化： .. automodule:: dspawpy.io.write :members: write_VESTA :noindex: volumetricData 指的是随空间位置变化的物理量，如电荷密度 rho，势能函数 potential，局域电荷密度 elf，部分电荷密度 pcharge，溶剂束缚电荷密度 rhoBound 等。这些数据在 DS-PAW 中以 volumetricData 类型保存。 - ``write_delta_rho_vesta`` 函数负责处理差分volumetricData可视化： .. automodule:: dspawpy.io.write :members: write_delta_rho_vesta :noindex: - ``average_along_axis`` 函数负责处理volumetricData面平均数据： .. automodule:: dspawpy.plot :members: average_along_axis :noindex: .. _ 能带数据处理: band能带数据处理 ============================================================================== .. admonition:: 知识点： 1. 脚本中调用get_band_data()读取数据，在读取数据时设置efermi=XX可将能量零点修改为指定值；设置zero_to_efermi=True, 可将能量零点修改问所读取的文件中的费米能级处 2. 脚本调用pymatgen的BSPlotter.get_plot()画图，在画图时可设置zero_to_efermi=True，将坐标轴能量零点修改到费米能级处。由于pymatgen在2023.8.17的一处关键更新将绘图函数返回对象从plt改成了axes，导致后续脚本可能出现不兼容。因此用户脚本相关部分已添加一个判断语句加以处理。 3. 能带两步算需从第一步自洽计算获取准确费米能级（从自洽获得的system.json），若获取失败，用户可在调用get_band_data读取数据时，利用参数efermi修改能量零点。例如：band_data=get_band_data('band.h5',efermi=-1.5) 4. 脚本中画图调用pymatgen中的BSPlotter.get_plot, 当判断体系为非金属时，设置zero_to_efermi会认为VBM为费米能级能量而非文件读取时的费米能级。故当体系为非金属时，在数据读取时设置zero_to_efermi=True和在画图时设置zero_to_efermi=True得到的图会有区别 执行本节所列的python脚本，程序会判断是否为金属体系。若为非金属体系，将要求选择是否平移费米能级到能量零点，请按提示操作。 普通能带处理 ----------------------------------------------------------- 参考 :download:`4bandplot.py <../dspawpy_proj/UserScripts/4bandplot.py>` ： .. literalinclude:: ../dspawpy_proj/UserScripts/4bandplot.py :language: python :linenos: .. admonition:: 知识点： 能带两步算需从自洽计算获取准确的费米能级（从system.json获取），若获取失败，用户可在 get_band_data 函数中指定 efermi 参数修改 执行代码可以得到类似以下能带图： .. figure:: ../imgs/us/4bandplot.png :height: 300px :align: center .. centered:: 二硫化钼能带示意图 将能带投影到每一种元素分别作图，数据点大小表示该元素对该轨道的贡献 ------------------------------------------------------------------------------------------- 参考 :download:`4bandplot_elt.py <../dspawpy_proj/UserScripts/4bandplot_elt.py>` ： .. literalinclude:: ../dspawpy_proj/UserScripts/4bandplot_elt.py :language: python :linenos: .. admonition:: 知识点： 1. 用户如果需要绘制能带投影的数据，此时需要使用 BSPlotterProjected模块 2. 使用 BSPlotterProjected模块中 get_elt_projected_plots 函数能够绘制每种元素对轨道贡献的能带图 执行代码可以得到类似以下能带图： .. figure:: ../imgs/us/4bandplot_elt.png :height: 300px :align: center .. centered:: 二硫化钼元素投影能带示意图 能带投影到不同元素的不同轨道 ---------------------------------------------------------- 参考 :download:`4bandplot_elt_orbit.py <../dspawpy_proj/UserScripts/4bandplot_elt_orbit.py>` ： .. literalinclude:: ../dspawpy_proj/UserScripts/4bandplot_elt_orbit.py :language: python :linenos: .. admonition:: 知识点： 1. 使用 BSPlotterProjected模块中 get_projected_plots_dots可以让用户来自定义需要绘制的某种元素某种轨道(L)的能带图 2. 例如 get_projected_plots_dots ({'Mo':['d'],'S':['s']})就是绘制Mo的d轨道和S的s轨道 执行代码可以得到类似以下能带图： .. figure:: ../imgs/us/4bandplot_elt_orbit.png :height: 500px :align: center .. centered:: 二硫化钼元素轨道投影能带示意图 将能带投影到不同原子的不同轨道 -------------------------------------------------------- 参考 :download:`4bandplot_patom_porbit.py <../dspawpy_proj/UserScripts/4bandplot_patom_porbit.py>` ： .. literalinclude:: ../dspawpy_proj/UserScripts/4bandplot_patom_porbit.py :language: python :linenos: .. admonition:: 知识点： 1. 使用 BSPlotterProjected模块中 get_projected_plots_dots_patom_pmorb 的自由度更高，可以让用户来自定义需要绘制的某些原子某些轨道的能带图 2. dictpa指定原子，dictio 指定该原子的轨道 3. 如果要将一些原子或一些轨道的投影分量叠加起来，请根据 get_projected_plots_dots_patom_pmorb 函数文档指定 sum_atoms 或 sum_morbs 参数 .. warning:: 1. 如果只选择单个轨道，且轨道名称不止一个字母（例如px、dxy、dxz），get_projected_plots_dots_patom_pmorb 函数将报错，详情见 `此处 `_ 执行代码可以得到类似以下能带图： .. figure:: ../imgs/us/4band_patom_porbit.png :height: 700px :align: center .. centered:: 二硫化钼原子投影能带示意图 能带反折叠处理 --------------------------------------------------------- 参考 :download:`4bandunfolding.py <../dspawpy_proj/UserScripts/4bandunfolding.py>` ： .. literalinclude:: ../dspawpy_proj/UserScripts/4bandunfolding.py :language: python :linenos: 执行代码可以得到类似以下能带图： .. figure:: ../imgs/us/4bandunfolding.png :height: 400px :align: center .. centered:: Cu3Au能带反折叠示意图 .. warning:: 此功能暂不支持将非金属材料的费米能级设为能量零点（默认价带顶为能量零点）。 band-compare能带对比图处理 --------------------------------------------------------- 将普通能带和wannier能带绘制在同一张图上 参考 :download:`4bandcompare.py <../dspawpy_proj/UserScripts/4bandcompare.py>` ： .. literalinclude:: ../dspawpy_proj/UserScripts/4bandcompare.py :language: python :linenos: 执行代码可得到类似如下能带对比曲线： .. figure:: ../imgs/phase4/wannier_band.png :height: 350px :align: center .. centered:: 晶体硅瓦尼尔能带与DFT能带示意图 .. dropdown:: :octicon:`code-review;1.5em;sd-text-light` API: get_band_data() :animate: fade-in-slide-down :color: info - ``get_band_data`` 函数负责读取能带数据如下： .. automodule:: dspawpy.io.read :members: get_band_data :noindex: .. _ 态密度数据处理: dos态密度数据处理 ============================================================================== 总的态密度 ------------------------------------------------------- 参考 :download:`5dosplot_total.py <../dspawpy_proj/UserScripts/5dosplot_total.py>` ： .. literalinclude:: ../dspawpy_proj/UserScripts/5dosplot_total.py :language: python :linenos: .. admonition:: 知识点： 1. 使用 get_dos_data 函数可以将DS-PAW计算得到的 dos.h5 文件转化为 pymatgen 支持的格式 2. 使用 DosPlotter模块获取到DS-PAW计算的 dos.h5 的数据 3. DosPlotter函数可以传递参数：stack参数表示画态密度是否加阴影， zero_at_efermi 表示是否在态密度图中进行将费米能量置零，这里设置 stack=False , zero_at_efermi=False 4. 使用 DosPlotter 模块中 add_dos 获取态密度的数据 5. DosPlotter模块中 get_plot函数 绘制态密度图 执行代码可以得到类似以下态密度图： .. figure:: ../imgs/us/5dos_total.png :height: 300px :align: center .. centered:: NiO态密度示意图 将态密度投影到不同的轨道上 -------------------------------------------------------- 参考 :download:`5dosplot_spd.py <../dspawpy_proj/UserScripts/5dosplot_spd.py>` ： .. literalinclude:: ../dspawpy_proj/UserScripts/5dosplot_spd.py :language: python :linenos: .. admonition:: 知识点： 使用 DosPlotter模块中 add_dos_dict 函数 获取投影态密度的数据，之后使用 get_spd_dos 将投影信息按照 spd 轨道投影输出 执行代码可以得到类似以下态密度图： .. figure:: ../imgs/us/5dos_spd.png :height: 300px :align: center .. centered:: NiO轨道投影态密度示意图 将态密度投影到不同的元素上 --------------------------------------------------------- 参考 :download:`5dosplot_elt.py <../dspawpy_proj/UserScripts/5dosplot_elt.py>` ： .. literalinclude:: ../dspawpy_proj/UserScripts/5dosplot_elt.py :language: python :linenos: .. admonition:: 知识点： 使用 DosPlotter模块中 add_dos_dict 函数 获取投影态密度的数据，之后使用 get_element_dos 将投影信息按照不同元素投影输出 执行代码可以得到类似以下态密度图： .. figure:: ../imgs/us/5dos_elt.png :height: 300px :align: center .. centered:: NiO元素投影态密度示意图 将态密度投影到不同原子的不同轨道上 --------------------------------------------------------- 参考 :download:`5dosplot_atom_orbit.py <../dspawpy_proj/UserScripts/5dosplot_atom_orbit.py>` ： .. literalinclude:: ../dspawpy_proj/UserScripts/5dosplot_atom_orbit.py :language: python :linenos: .. admonition:: 知识点： 1. 使用 get_site_orbital_dos函数提取dos数据中特定原子特定轨道的贡献，dos_data.structure[0],Orbital(4) 表示获取第1个原子的dxy轨道的态密度，get_site_orbital_dos函数中序号从0开始 2. 运行此脚本，根据提示选择元素和轨道，可以得到相应的态密度图 执行代码可以得到类似以下态密度图： .. figure:: ../imgs/us/5dos_atom_orbit.png :height: 300px :align: center .. centered:: NiO原子轨道态密度示意图 将态密度投影到不同原子的分裂d轨道(t2g, eg)上 -------------------------------------------------------- 参考 :download:`5dosplot_t2g_eg.py <../dspawpy_proj/UserScripts/5dosplot_t2g_eg.py>` ： .. literalinclude:: ../dspawpy_proj/UserScripts/5dosplot_t2g_eg.py :language: python :linenos: .. admonition:: 知识点： 1. 使用 get_site_t2g_eg_resolved_dos函数 提取dos数据中特定原子的 t2g和 eg轨道的贡献，这是获取第2个原子的t2g和eg轨道的贡献 2. 运行此脚本，根据提示选择原子序号，可以得到相应的态密度图 执行代码可以得到类似以下态密度图： .. figure:: ../imgs/us/5dos_t2g_eg.png :height: 300px :align: center .. centered:: NiO分裂d轨道原子投影态密度示意图 .. admonition:: 知识点： 如果元素不含d轨道，会画成空白图片 d-带中心分析 --------------------------------------------------------- 以 Pb-slab 体系为例，对 Pt 原子进行 d-band 中心分析： 参考 :download:`5center_dband.py <../dspawpy_proj/UserScripts/5center_dband.py>` ： .. literalinclude:: ../dspawpy_proj/UserScripts/5center_dband.py :language: python :linenos: 执行代码可以得到类似以下结果： .. code-block:: python spin=1 -1.785319344084034 .. note:: 目前仅支持全部原子平均的 d 轨道中心，不支持元素、原子投影或其他轨道，也不支持选择自旋方向或能量范围。 ``get_dos_data`` 函数负责处理态密度数据： .. dropdown:: :octicon:`code-review;1.5em;sd-text-light` API: get_dos_data() :animate: fade-in-slide-down :color: info .. automodule:: dspawpy.io.read :members: get_dos_data :noindex: .. _ 能带和态密度共同显示: bandDos能带和态密度共同显示 ============================================================================== 以应用教程中Si体系为例： 将能带和态密度显示在一张图上 ------------------------------------------------------------------------------ 参考 :download:`6bandDosplot.py <../dspawpy_proj/UserScripts/6bandDosplot.py>` ： .. literalinclude:: ../dspawpy_proj/UserScripts/6bandDosplot.py :language: python :linenos: 执行代码可以得到类似以下能带态密度图： .. figure:: ../imgs/us/6bandDos.png :height: 500px :align: center .. centered:: 单晶硅能带-态密度示意图 将能带和投影态密度显示在一张图上 ------------------------------------------------------------------------------ 参考 :download:`6bandPdosplot.py <../dspawpy_proj/UserScripts/6bandPdosplot.py>` ： .. literalinclude:: ../dspawpy_proj/UserScripts/6bandPdosplot.py :language: python :linenos: 执行代码可以得到类似以下能带态密度图： .. figure:: ../imgs/us/6bandPdos.png :height: 500px :align: center .. centered:: 单晶硅能带-投影态密度示意图 .. warning:: 1. 给定投影能带数据，默认将沿着元素投影；给定普通能带数据（或体系所含元素种类超过4种），则不投影，并输出警告 2. 给定投影态密度数据，默认也将沿着元素投影，可以换成沿着轨道投影，或者不投影；给定普通态密度数据且未关闭态密度投影选项 BSDOSPlotter(dos_projection=None)，pymatgen绘图程序将报错，这就是上面特意准备了一个 6bandDosplot.py 文件的原因。 .. dspawpy还提供了一个优化后的band-compare能带对比图绘图函数 ``pltbd`` .. .. dropdown:: :octicon:`code-review;1.5em;sd-text-light` API: pltbd() .. :animate: fade-in-slide-down .. :color: info .. .. automodule:: dspawpy.plot .. :members: pltbd .. :noindex: .. 新老函数绘图效果对比： .. .. figure:: ../imgs/newbd.png .. :height: 500px .. :align: center .. .. centered:: 新 pltbd 绘图效果 .. .. figure:: ../imgs/oldbd.png .. :height: 500px .. :align: center .. .. centered:: 旧 get_plot 绘图效果 .. _ 光学性质数据处理: optical光学性质数据处理 ============================================================================== 以快速入门Si体系光学性质计算得到的 :guilabel:`scf.h5` 为例（注：输出文件名与task一致，task = scf，io.optical = true可计算光学性质）： 反射率数据处理，参考 :download:`7optical.py <../dspawpy_proj/UserScripts/7optical.py>` ： .. literalinclude:: ../dspawpy_proj/UserScripts/7optical.py :language: python :linenos: .. admonition:: 知识点： Reflectance为光学性质中的一种，用户可以根据自己的需求将该关键词修改为“AbsorptionCoefficient”或”ExtinctionCoefficient”或”RefractiveIndex”，分别对应吸收系数、消光系数和折射率 执行代码可以得到类似以下反射率随能量变化的曲线： .. figure:: ../imgs/us/7optical.png :height: 300px :align: center .. centered:: 单晶硅Reflectance随光子能量变化趋势示意图 .. dropdown:: :octicon:`code-review;1.5em;sd-text-light` API: plot_optical() :animate: fade-in-slide-down :color: info .. automodule:: dspawpy.plot :members: plot_optical :noindex: .. _ 过渡态计算数据处理: neb过渡态计算数据处理 ============================================================================== 以快速入门 H 在 Pt(100) 表面扩散为例进行介绍： 输入文件之生成中间构型 ---------------------------------------------------------- - 参考 :download:`8neb_interpolate_structures.py <../dspawpy_proj/UserScripts/8neb_interpolate_structures.py>` ： .. literalinclude:: ../dspawpy_proj/UserScripts/8neb_interpolate_structures.py :language: python :linenos: .. admonition:: 知识点： 1. 用户可以根据需要自行修改插点个数，设置为8得到包含8个结构文件的文件夹，其中中间构型为6个 2. neb.linear_interpolate为线性插值方法，若使用idpp插值请注释该行放开下一行 绘制能垒图 ---------------------------------------------------------- neb.iniFin = true/false ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 当设置 neb.iniFin = true/false 时，都可使用读取neb计算的路径进行势垒分析（确保初末态计算文件在neb计算路径下）： - 参考 :download:`8neb_barrier_CubicSpline.py <../dspawpy_proj/UserScripts/8neb_barrier_CubicSpline_directory.py>` ： .. literalinclude:: ../dspawpy_proj/UserScripts/8neb_barrier_CubicSpline_directory.py :language: python :linenos: .. note:: 运行上述脚本后，可得到类似以下三次样条插值的势垒曲线： .. figure:: ./../imgs/phase2/neb-barrier_cs.png :align: center :height: 500px 对于这个算例，三次样条插值后，曲线会出现不合理的“下坠”区域，这是三次样条插值算法的特性所决定的。 dspawpy 内部调用 `scipy 的插值算法 `_ ，上面这个脚本我们以三次样条插值算法为例，它在 scipy 文档中的定义为： .. code-block:: python class scipy.interpolate.CubicSpline(x, y, axis=0, bc_type='not-a-knot', extrapolate=None) 关键词参数包括 axis, bc_type, extrapolate，具体含义见 `scipy.interpolate.CubicSpline `_ 。我们可以往 plot_barrier 函数中指定相应的关键词参数（axis, bc_type, extrapolate），将其传递给 scipy.interpolate.CubicSpline 这个类处理。 下面我们采用 :download:`8neb_barrier.py <../dspawpy_proj/UserScripts/8neb_barrier.py>` 脚本，对比三种算法插值绘制出的曲线： .. literalinclude:: ../dspawpy_proj/UserScripts/8neb_barrier.py :language: python :linenos: :emphasize-lines: 15, 21, 23 .. figure:: ./../imgs/phase2/neb-barrier.png :height: 500px :align: center .. admonition:: 知识点： 1. 选择合适的插值算法对于优化最终曲线的呈现效果至关重要 2. 多数情况下，选择 pchip 这种单调三次样条插值算法即可达到较好效果，这也是默认调用的插值算法 neb.iniFin = true ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 当设置 neb.iniFin = true 时，读取neb计算所得的 neb.h5/neb.json 文件可快速进行势垒分析： - 参考 :download:`8neb_barrier_CubicSpline.py <../dspawpy_proj/UserScripts/8neb_barrier_CubicSpline.py>` ： .. literalinclude:: ../dspawpy_proj/UserScripts/8neb_barrier_CubicSpline.py :language: python :linenos: :emphasize-lines: 6 处理得到的势垒图与之前读取路径效果一致。 .. note:: 1. neb.h5 和 neb.json文件所存能量为TotalEnergy，如需得精确的势垒值，建议采用读取neb计算路径的方法进行处理（取TotalEnergy0） 过渡态计算数据处理 ---------------------------------------------------------- NEB计算完成后，一般要画出能垒图观察，并检查各插值构型最终的受力大小，从而确保受力小于指定的阈值。如果结果异常，还需要检查各插值构型在结构优化过程中的受力与能量变化趋势，判断是否真正“收敛”。上述这些操作至少需要三次循环，为简化操作，我们提供了一个一步到位的总结函数 ``summary`` ： - 参考 :download:`8neb_check_results.py <../dspawpy_proj/UserScripts/8neb_check_results.py>` ： .. literalinclude:: ../dspawpy_proj/UserScripts/8neb_check_results.py :language: python :linenos: .. admonition:: 知识点： 1. 此脚本将以表格形式打印各构型的能量和受力、绘制能垒图，并绘制中间构型的能量与受力收敛过程图 2. 若 ``neb.iniFin = false`` ，用户必须将自洽计算的结果文件 :guilabel:`scf.h5` 或 :guilabel:`system.json` 文件复制到对应的初末态子文件夹中，否则程序无法读取初末态能量和受力信息，将报错退出 3. 默认情况下，能垒图存储在neb计算目录外层，各中间构型的能量与受力收敛图存储在各子文件夹 执行代码可以得到类似如下NEB计算各构型的能量和受力表格： +-----+-------------+-----------------------+---------------------+---------------------------+ |Image|Force (eV/Å) |Reaction coordinate (Å)| Energy (eV) | Delta energy (eV) | +-----+-------------+-----------------------+---------------------+---------------------------+ |00 | 0.1803 | 0.0000 | -39637.0984 | 0.0000 | +-----+-------------+-----------------------+---------------------+---------------------------+ |01 | 0.0263 | 0.5428 | -39637.0186 | 0.0798 | +-----+-------------+-----------------------+---------------------+---------------------------+ |02 | 0.0248 | 1.0868 | -39636.8801 | 0.2183 | +-----+-------------+-----------------------+---------------------+---------------------------+ |03 | 0.2344 | 1.5884 | -39636.9984 | 0.1000 | +-----+-------------+-----------------------+---------------------+---------------------------+ |04 | 0.0141 | 2.0892 | -39637.0900 | 0.0084 | +-----+-------------+-----------------------+---------------------+---------------------------+ 除了可以获得能垒图外，还可以得到各中间构型的能量和受力收敛曲线（以02构型为例） .. figure:: ./../imgs/phase2/neb-energy.png :height: 500px :align: center .. _观察NEB链: 观察NEB链 ---------------------------------------------------------- 此处所说NEB链指的是各插值构型（structure00.as, structure01.as, ...）之间的几何关系，而不是单个构型在优化过程中的变化。 - NEB任务计算量较大，观察NEB链有助于判断NEB计算的收敛速度。另外，通过插值算法生成中间结构后，经常需要预览NEB链。这些需求可以通过 :download:`8neb_visualize.py <../dspawpy_proj/UserScripts/8neb_visualize.py>` 脚本实现： .. literalinclude:: ../dspawpy_proj/UserScripts/8neb_visualize.py :language: python :linenos: .. admonition:: 知识点： 1. 此脚本生成neb_movie*.json文件后，通过 ``Device Studio`` --> ``Simulator`` --> ``DS-PAW`` --> ``Analysis Plot`` 打开 json 文件即可观察 2. directory指定为NEB计算主路径，需提供neb计算完成后的完整文件夹 3. 该脚本支持处理正在进行中的（即未完成的）neb计算文件，方便用户监测实时轨迹 4. xyz文件可通过 OVITO 软件打开查看：通过 ``Device Studio`` --> ``Simulator`` --> ``OVITO`` 打开可视化界面，将xyz文件拖入即可 5. 结构信息读取优先级：latestStructureXX.as > h5 > json；设置ignorels为True后，先尝试读取 h5 中的数据，失败则读取 json 中的数据 计算构型间距 ---------------------------------------------------------- - 参考这个脚本 :download:`8calc_dist.py <../dspawpy_proj/UserScripts/8calc_dist.py>` ： .. literalinclude:: ../dspawpy_proj/UserScripts/8calc_dist.py :language: python :linenos: neb续算 -------------------------------------------------------- - 如果需对neb进行续算，可参考 :download:`8neb_restart.py <../dspawpy_proj/UserScripts/8neb_restart.py>` ： .. literalinclude:: ../dspawpy_proj/UserScripts/8neb_restart.py :language: python :linenos: 具体效果见 :ref:`neb过渡态计算续算说明`。 .. dropdown:: :octicon:`code-review;1.5em;sd-text-light` API: write_neb_structures(), plot_barrier(), summary(), get_distance(), write_movie_json(), write_xyz(), restart() :animate: fade-in-slide-down :color: info - ``write_neb_structures`` 函数负责生成中间构型： .. automodule:: dspawpy.diffusion.neb :members: write_neb_structures :noindex: - ``plot_barrier`` 函数负责绘制能垒图： .. automodule:: dspawpy.diffusion.nebtools :members: plot_barrier :noindex: - ``summary`` 函数负责总结NEB计算任务的说明文档： .. automodule:: dspawpy.diffusion.nebtools :members: summary :noindex: - ``get_distance`` 函数可以计算两个构型间的距离： .. automodule:: dspawpy.diffusion.nebtools :members: get_distance :noindex: - ``write_movie_json`` 和 ``write_xyz`` 函数可以将中间构型写入json或者xyz文件： .. automodule:: dspawpy.diffusion.nebtools :members: write_movie_json, write_xyz :noindex: - ``restart`` 函数负责重启NEB计算： .. automodule:: dspawpy.diffusion.nebtools :members: restart :noindex: .. _ 声子计算数据处理: phonon声子计算数据处理 ============================================================================== 以MgO体系的声子能带态密度计算得到的 :guilabel:`phonon.h5` 为例： 如果未安装 phonopy，运行下列脚本时会弹出 ``no module named 'phonopy'`` 信息，不影响程序正常运行 声子能带数据处理 ----------------------------------------------------- - 参考 :download:`9phonon_bandplot.py <../dspawpy_proj/UserScripts/9phonon_bandplot.py>` ： .. literalinclude:: ../dspawpy_proj/UserScripts/9phonon_bandplot.py :language: python :linenos: 执行代码可以得到类似以下声子能带曲线： .. figure:: ../imgs/phase3/phonon-nonac.png :height: 500px :align: center 声子态密度数据处理 ---------------------------------------------------- - 参考 :download:`9phonon_dosplot.py <../dspawpy_proj/UserScripts/9phonon_dosplot.py>` ： .. literalinclude:: ../dspawpy_proj/UserScripts/9phonon_dosplot.py :language: python :linenos: 执行代码可以得到类似以下声子态密度曲线： .. figure:: ../imgs/us/9phonon_dosplot.png :height: 350px :align: center .. centered:: 单晶硅声子态密度示意图 声子热力学数据处理 ----------------------------------------------------- 可以参考 :download:`9phonon_thermal.py <../dspawpy_proj/UserScripts/9phonon_thermal.py>` ： .. literalinclude:: ../dspawpy_proj/UserScripts/9phonon_thermal.py :language: python :linenos: 执行代码可以得到类似以下声子热力学曲线： .. figure:: ../imgs/us/9phonon.png :height: 500px :align: center .. centered:: 单晶硅声子热力学性质示意图 .. dropdown:: :octicon:`code-review;1.5em;sd-text-light` API: get_phonon_band_data(), get_phonon_dos_data(), plot_phonon_thermal() :animate: fade-in-slide-down :color: info - ``get_phonon_band_data`` 函数负责读取声子能带： .. automodule:: dspawpy.io.read :members: get_phonon_band_data :noindex: - ``get_phonon_dos_data`` 函数负责读取声子态密度： .. automodule:: dspawpy.io.read :members: get_phonon_dos_data :noindex: - ``plot_phonon_thermal`` 函数负责绘制声子热力学性质图： .. automodule:: dspawpy.plot :members: plot_phonon_thermal :noindex: .. _ 分子动力学模拟数据处理: aimd分子动力学模拟数据处理 ============================================================================== 以快速入门 :math:`H_{2}O` 分子体系分子动力学模拟得到的 :guilabel:`aimd.h5` 为例： 轨迹文件转换格式为.xyz或.dump ------------------------------------------------------ 从aimd输出的hdf5文件中读取数据，并生成轨迹文件 生成的.xyz或.dump格式文件，可拖入OVITO观察，通过 Device Studio --> Simulator --> OVITO 打开 OVITO 可视化界面，将xyz文件或dump文件拖入OVITO即可。 参考 :download:`10write_aimd_traj.py <../dspawpy_proj/UserScripts/10write_aimd_traj.py>` ： .. literalinclude:: ../dspawpy_proj/UserScripts/10write_aimd_traj.py :language: python :linenos: 执行代码将生成.xyz和.dump格式的轨迹文件，该文件可通过OVITO打开。关于结构文件转化的更多细节，请参考 :ref:`结构转化` .. admonition:: 知识点： OVITO 与 dspawpy 都不支持将非正交晶胞的体系保存为dump文件 动力学过程中能量、温度等变化曲线 ---------------------------------------------------- - 参考 :download:`10check_aimd_conv.py <../dspawpy_proj/UserScripts/10check_aimd_conv.py>` ： .. literalinclude:: ../dspawpy_proj/UserScripts/10check_aimd_conv.py :language: python :linenos: 执行代码将生成如下组合图： .. figure:: ../imgs/phase2/aimd-joined.png :height: 800px :align: center 均方位移（MSD） 分析 ------------------------------------------------------ - 参考 :download:`10aimd_msd.py <../dspawpy_proj/UserScripts/10aimd_msd.py>` ： .. literalinclude:: ../dspawpy_proj/UserScripts/10aimd_msd.py :language: python :linenos: 执行代码将生成类似如下图片： .. figure:: ../imgs/us/10MSD.png :height: 500px :align: center .. centered:: 均方位移（MSD）示意图 均方根偏差（RMSD） 分析 ------------------------------------------------------ - 参考 :download:`10aimd_rmsd.py <../dspawpy_proj/UserScripts/10aimd_rmsd.py>` ： .. literalinclude:: ../dspawpy_proj/UserScripts/10aimd_rmsd.py :language: python :linenos: 执行代码将生成类似如下图片： .. figure:: ../imgs/us/10RMSD.png :height: 500px :align: center .. centered:: 均方根偏差（RMSD）示意图 原子对分布函数或径向分布函数（RDFs） 分析 ------------------------------------------------------ - 参考 :download:`10aimd_rdf.py <../dspawpy_proj/UserScripts/10aimd_rdf.py>` ： .. literalinclude:: ../dspawpy_proj/UserScripts/10aimd_rdf.py :language: python :linenos: 执行代码将生成类似如下图片： .. figure:: ../imgs/us/10RDF.png :height: 500px :align: center .. centered:: 径向分布函数（RDFs）示意图 - 这部分涉及的统计学计算较复杂，更多细节请参考函数API .. dropdown:: :octicon:`code-review;1.5em;sd-text-light` API: plot_aimd(), get_lagtime_msd(), plot_msd(), get_rs_rdfs(), plot_rdf(), get_lagtime_rmsd(), plot_rmsd() :animate: fade-in-slide-down :color: info - ``plot_aimd`` 函数可用于协助检查AIMD计算过程中关键物理量的收敛过程： .. automodule:: dspawpy.plot :members: plot_aimd :noindex: - ``get_*`` 和 ``plot_*`` 函数负责读取AIMD计算过程中的关键物理量： .. automodule:: dspawpy.analysis.aimdtools :members: get_lagtime_msd, plot_msd, get_rs_rdfs, plot_rdf, get_lagtime_rmsd, plot_rmsd :noindex: - 提取数据自行处理请参考： .. code-block:: python :linenos: from dspawpy.io.read import get_sinfo from dspawpy.io.structure import read aimd_h5_files = ['aimd1.h5','aimd2.h5','aimd3.h5'] # 可以从计算完成的多个 aimd.h5 中依次抽取数据并合并 # 一次性读取多个aimd.h5文件中的数据并创建pymatgen的Structures列表 pymatgen_structures = read(datafile=aimd_h5_files) # 或者将数组提取出来 for i, df in enumerate(aimd_h5_files): # 依次获取每个aimd.h5文件的数据 Nstep, elements, positions, lattices, D_mag_fix = get_sinfo(df) # Nstep 表示总离子步数 （int） # elements 表示元素列表，（Natom x 1） # positions 表示离子坐标，（Nstep x Natom x 3） # lattices 表示晶胞矩阵，（Nstep x 3 x 3 # D_mag_fix 磁矩、自由度相关信息字典 .. _ Polarization铁电极化数据处理: Polarization铁电极化数据处理 ============================================================================== 以快速入门 :math:`HfO_{2}` 体系铁电计算得到的系列 :guilabel:`scf.h5` 为例： - 参考 :download:`11Ferri.py <../dspawpy_proj/UserScripts/11Ferri.py>` ： .. literalinclude:: ../dspawpy_proj/UserScripts/11Ferri.py :language: python :linenos: 执行代码将生成如下组合图： .. figure:: ../imgs/us/11pol.png :height: 500px :align: center .. centered:: 12组结构对应极化数值 查看首尾构型的铁电极化数值，可以参考如下： .. code-block:: python :linenos: from dspawpy.plot import plot_polarization_figure plot_polarization_figure(directory='.', annotation=True, annotation_style=1) # 显示首尾构型的铁电极化数值 执行代码将生成如下组合图： .. figure:: ../imgs/phase3/Ferri-Pola-anno1.png :height: 500px :align: center .. centered:: 12组结构对应极化数值（带首尾构型数值） 也可以使用第二种批注样式： .. code-block:: python :linenos: from dspawpy.plot import plot_polarization_figure plot_polarization_figure(directory='.', annotation=True, annotation_style=2) # 显示首尾构型的铁电极化数值 执行代码将生成如下组合图： .. figure:: ../imgs/phase3/Ferri-Pola-anno2.png :height: 500px :align: center .. centered:: 12组结构对应极化数值（带首尾构型数值） .. dropdown:: :octicon:`code-review;1.5em;sd-text-light` API: plot_polarization_figure() :animate: fade-in-slide-down :color: info - ``plot_polarization_figure`` 函数负责绘制铁电极化图： .. automodule:: dspawpy.plot :members: plot_polarization_figure :noindex: .. _ ZPE零点振动能数据处理: ZPE零点振动能数据处理 ============================================================================== 以快速入门 :math:`CO` 体系频率计算得到的 :guilabel:`frequency.txt` 文件为例，计算零点振动能，基于以下公式： .. math:: ZPE=\sum_{i=1}^{3 N} \frac{h \nu_i}{2} 其中， :math:`\nu_i` 是简正频率， :math:`h` 是普朗克常数（ :math:`6.626\times10^{-34} J\cdot s`）， :math:`N` 是原子数。 - 参考 :download:`12getZPE.py <../dspawpy_proj/UserScripts/12getZPE.py>` ： .. literalinclude:: ../dspawpy_proj/UserScripts/12getZPE.py :language: python :linenos: 执行代码结果文件将保存到 :guilabel:`ZPE.dat` 文件中，文件内容如下： .. code-block:: text Data read from D:\quickstart\2.13\frequency.txt Frequency (meV) 284.840038 --> Zero-point energy, ZPE (eV): 0.142420019 .. dropdown:: :octicon:`code-review;1.5em;sd-text-light` API: getZPE() :animate: fade-in-slide-down :color: info - ``getZPE`` 函数负责计算零点振动能： .. automodule:: dspawpy.io.utils :members: getZPE :noindex: .. _ TS的热校正能: TS的热校正能 ============================================================================== 吸附质的熵变对能量的贡献 ----------------------------------------------------- 计算基于以下公式： .. math:: \Delta S_{a d s}\left(0 \rightarrow T, P^0\right)=S_{v i b}=\sum_{i=1}^{3 N}\left[\frac{N_{\mathrm{A}} h \nu_i}{T\left(e^{h \nu_i / k_{\mathrm{B}} T}-1\right)}-R \ln \left(1-e^{-h \nu_i / k_{\mathrm{B}} T}\right)\right] 其中， :math:`\Delta S_{a d s}` 表示吸附质的熵变，根据简谐近似计算。:math:`S_{v i b}` 表示振动熵。 :math:`\nu_i` 是简正频率，:math:`N_A` 是阿伏伽德罗常数（ :math:`6.022\times10^{23} mol^{-1}` ）， :math:`h` 是普朗克常数（ :math:`6.626\times10^{-34} J\cdot s`）， :math:`k_B` 是玻尔兹曼常数（ :math:`1.38\times10^{-23} J\cdot K^{-1}`）， :math:`R` 是理想气体常数（ :math:`8.314 J\cdot mol^{-1}\cdot K^{-1}`）， :math:`T` 是体系温度，单位 :math:`K`。 - 参考 :download:`13getTSads.py <../dspawpy_proj/UserScripts/13getTSads.py>` ： .. literalinclude:: ../dspawpy_proj/UserScripts/13getTSads.py :language: python :linenos: 执行代码结果文件将保存到 :guilabel:`TS.dat` 文件中，文件内容如下： .. code-block:: text Data read from D:\quickstart\2.13\frequency.txt Frequency (THz) 68.873994 --> Entropy contribution, T*S (eV): 4.7566990201851275e-06 理想气体的熵变对能量的贡献 ----------------------------------------------------- 计算基于如下公式： .. math:: \begin{aligned} S(T, P) & =S\left(T, P^{\circ}\right)-k_{\mathrm{B}} \ln \frac{P}{P^{\circ}} \\ & =S_{\text {trans }}+S_{\text {rot }}+S_{\text {elec }}+S_{\text {vib }}-k_{\mathrm{B}} \ln \frac{P}{P^{\circ}} \end{aligned} 其中： .. math:: S_{\text {trans }}=k_{\mathrm{B}}\left\{\ln \left[\left(\frac{2 \pi M k_{\mathrm{B}} T}{h^2}\right)^{3 / 2} \frac{k_{\mathrm{B}} T}{P^{\circ}}\right]+\frac{5}{2}\right\} .. math:: S_{\mathrm{rot}}= \begin{cases}0 & \text {, monatomic } \\ k_{\mathrm{B}}\left[\ln \left(\frac{8 \pi^2 I k_{\mathrm{B}} T}{\sigma h^2}\right)+1\right] & \text {, linear } \\ k_{\mathrm{B}}\left\{\ln \left[\frac{\sqrt{\pi I_A I_{\mathrm{B}} I_{\mathrm{C}}}}{\sigma}\left(\frac{8 \pi^2 k_B T}{h^2}\right)^{3 / 2}\right]+\frac{3}{2}\right\} & \text {, nonlinear }\end{cases} .. math:: S_{\text {elec }}=k_{\mathrm{B}} \ln [2 \times(\text { total spin })+1] .. math:: S_{\text {vib }}=k_{\mathrm{B}} \sum_i^{\text {vib DOF }}\left[\frac{\epsilon_i}{k_{\mathrm{B}} T\left(e^{\epsilon_i / k_{\mathrm{B}} T}-1\right)}-\ln \left(1-e^{-\epsilon_i / k_{\mathrm{B}} T}\right)\right] .. math:: S_{\text {elec }}=k_{\mathrm{B}} \ln [2 \times(\text { total spin })+1] 其中：:math:`I_A` 到 :math:`I_C` 是非线性分子的三个主惯性矩， :math:`I` 是线性分子的简并惯性矩， :math:`\sigma` 是分子的对称数。另外，momatomic 表示单原子分子，linear 表示线性分子，nonlinear 表示非线性分子。total spin 是总自旋数。 vib DOF 表示振动自由度。 - 参考 :download:`13getTSgas.py <../dspawpy_proj/UserScripts/13getTSgas.py>` 脚本处理： .. literalinclude:: ../dspawpy_proj/UserScripts/13getTSgas.py :language: python :linenos: .. dropdown:: :octicon:`code-review;1.5em;sd-text-light` API: getTSads(), getTSgas() :animate: fade-in-slide-down :color: info - ``getTSads`` 函数负责计算吸附质熵变对能量的贡献： .. automodule:: dspawpy.io.utils :members: getTSads :noindex: - ``getTSgas`` 函数负责计算理想气体熵变对能量的贡献： .. automodule:: dspawpy.io.utils :members: getTSgas :noindex: 附录 ============================================================================== - 快速下载所有脚本，点击 :download:`用户脚本.zip <../dspawpy_proj/UserScripts.zip>` - `dspawpy更新日志 `_ .. - 更多dspawpy的细节请移步 `API `_